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^ ! Abstract 

We introduced a method to obtain the continuum description of the elastic properties of mono- 



id 



layer h-BN through ab initio density functional theory. This thermodynamically rigorous contin- 
uum description of the elastic response is formulated by expanding the elastic strain energy density 
in a Taylor series in strain truncated after the fifth-order term, we obtained a total of fourteen 
nonzero independent elastic constants for the up to tenth-order tensor. We predicted the pressure 
dependent second-order elastic moduli. This continuum formulation is suitable for incorporation 

into the finite element method. 

PACS numbers: 62.25.-g, 81.40.Jj, 71.15.Mb, 71.15.Nc 
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I. INTRODUCTION 



The area of research on 2D nanomaterials with potential next generation device appli- 
cation has seen tremendous progress in the past few recent years. An example of such 
nanostructures is hexagonal boron nitride (h-BN) monolayer which is analog of graphene 
having a honeycomb lattice structure [l|. Hexagonal boron nitride is chemically inert, has 
a high thermal conductivity, and is highly temperature resistant to oxidation. Due to its 
outstanding properties, h-BN has found wide applications in micro and nano-devices such as 
insulator with high thermal conductivity in electronic devices |l|, ultraviolet-light emitter in 
optoelectronics 2j-|6j, and as nano- fillers in high strength and thermal conductive nanocom- 
posites [3, [8j . In very recent works, it has been shown that a tunable band gap nanosheet 
can be constructed by fabrication of hybrid nanostructures made of graphene/h-BN domains 
which opens a new venue for huge research in the application of h-BN for electronics [9|, Il0| . 

Given the aforementioned potential applications, the complete knowledge of mechanical 
and physical properties of h-BN monolayer, however, is still lacking. The previous primary 
works have reported that h-BN monolayer has a bulk modulus around 160 Pa.m and a 



bending modulus around 31.2GPa (Ref. [l|, Illl-ll3|) (- Amir, could you please double check 
this? It does not make sense.) Several experimental and atomistic simulation studies, mostly 
on graphene and carbon nanotubes, have probed that 2D nanosheets and nanotubes usually 
show a nonlinear elastic deformation during the tension up to the intrinsic strength of the 
material followed by a strain softening up to the fracture 14-ll8|. To establish a continuum 
based framework to capture this nonlinear elastic behavior of the 2D nanosheets, the hi ghe r 
order of the elastic constants must be considered in the strain energy density function 19ll20|. 
In such a model, the strain energy density is expanded in a Taylor series to include both 
quadratic as well as higher order terms in strain. The quadratic term accounts for the 
linear elastic response of the material while the cubic and higher order terms account the 
strain softening of the elastic stiffness. The higher order terms also can be used to define 
other anharmonic properties of this 2D nanostructure including phenomena such as thermal 
expansion, phonon-phonon interaction, etc 21 ] . 

The goal of this paper is to find the continuum description of the elastic properties of 
monolayer h-BN. To achieve that, we first examine the elastic properties of h-BN monolayers 
using ab initio density functional theory. We adopt a fifth-order series expansion of the 
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FIG. 1. Atomic super cell (24 atoms) of h-BN in the undefbrmed reference configuration. 



strain energy density function in order to model the inplane elastic properties of h-BN 
and demonstrate that the resulting continuum description now with fourteen independent 
elastic constants describes accompanying ab initio DFT calculations with high accuracy in 
the infinitesimal strain regime as well as at finite strains, including the strain at the intrinsic 
stress and beyond. A higher rank tensor is associated with each term of the series expansion 
and the components of the tensor represent the continuum elastic properties. Previous 
authors had determined the nonzero independent tensor components that correspond to the 
symmetry elements of graphene for the second-, third-, fourth-order terms and fifth-order 
term from stress-strain response for graphene 20J. We extended the method with least- 
squares solution to over-determined (up to eighth-rank tensor) and well-determined (tenth- 
rank tensor) linear equations. We applied this advanced method to obtain the continuum 
description of the elastic properties of of monolayer h-BN in the following sections. 



II. NONLINEAR ELASTICITY THEORY APPLIED TO 2D HEXAGONAL STRUC- 
TURE 



We used a super cell containing 12 B and 12 N atoms in one plane, with periodic boundary 
conditions. The undeformed reference configuration is shown in Fig. [H with lattice vectors 
Hi, i = 1,2,3. When a macroscopically homogeneous deformation (deformation gradient 
tensor 22J F) applied, the lattice vectors of the deformed h-BN are hj = FHj. The La- 



grangian strain [23] is defined as T7=-(F T F — I), where I is the identity tensor. The strain 
energy density has functional form of $ = &(tj)- The elastic properties of a material are 
determined from $, which is quadratic in strain for a linear elastic material. Nonlinear elas- 



tic constitutive behavior is established by expanding $ in a Taylor series in terms of powers 
of strain 77. The symmetric second Piola-Kirchhoff stress tensor, Ey, can be expressed (up 



to fifth order) as [20|: 
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where rjij is Lagrangian elastic strain. Summation convention is employed for repeating 
indices; lower case subscripts range from 1 to 3. Herein C denotes each higher-order elastic 
modulus tensor; the rank of each tensor corresponds to the number of subscripts. The 
second-order elastic constants (SOEC), Cyjw, third-order elastic constants (TOEC), Cijki mn , 
fourth-order elastic constants (FOEC), Cijkimnop, and fifth- order elastic constants (FFOEC), 
Cijkimnopqr, are given by the components of the fourth-, sixth-, eighth-, and tenth-rank 
tensors, respectively. 



We used conventional Voigt notation (24j for subscripts: 11 — > 1, 22 — > 2, 33 — > 3, 23 — > 
4, 31 — > 5, and 12 — > 6. Please note that for strain 774 = 27723, 775 = 2r/ 31 , rje = 277 12 . Eqs. (1) 
can be rewritten as 
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+ ^CijklVjVkVl + -tjCuklmVjVkVlVm- (2) 

where the summation convention for upper case subscripts runs from 1 to 6. 

In this study, we modeled the monolayer h-BN as two dimensional (2D) structure and 
assume that the deformed state of the monolayer h-BN is such that the contribution of bend- 
ing to the strain energy density is negligible as compared to the in-plane strain contribution. 
This assumption is reasonable since the radius of curvature of out-of-plane deformation are 
significantly larger than the in-plane inter-atomic distance. Then the stress state of mono- 
layer h-BN under those assumptions can be assumed to be 2D and we only consider the 
in-plane stress and strain components for these kind of structures. 

The components of the TOEC, FOEC, and FFOEC tensors can be determined based on 



the symmetries of the graphene atomic lattice point group Dgh which con- sists of a sixfold 
rotational axis and six mirror planes as formulated in ref [20j. 

The fourteen independent elastic constants of h-BN are determined by a least-squares 
fit to stress-strain results from ab initio DFT simulations in two steps. At the first step, 
we use least-squares fit to five stress-strain responses. Five relationships between stress and 
strain are necessary because there are five independent FFOECs. We obtain the stress-strain 
relationships by simu- lating the following deformation states: uniaxial strain in the zigzag 
direction; uniaxial strain in the armchair direction; and, equibiaxial strain. From the first 
step, we the components of SOEC, TOEC, FOEC are over-determined (i.e, the number of 
linearly independent variables are greater than the number of constrains), and the FFOEC 
are well-determined (the number of linearly independent variables are equal to the number 
of constrains). Under such circumstance, the second step is needed: least-square solution to 
these over- and well- determined linear equations. 

At the first step, we carried out three deformations, uniaxial strain in the zigzag direction 
(case z), uniaxial strain in the armchair direction (case a) and equibiaxial strain (case b). 
For uniaxial strain in the zigzag direction, the strain tensor is, 



v!j = o r] z o 


where rj z is the amount of strain in zigzag direction. 

For a given strain tensor, the associated deformation gradient tensor is not unique. The 
various possible solutions differing from one to another by a rigid rotation. Here the lack of 
a one-to-one map relationship between the strain tensor and deformation gradient tensor is 
not concern since the calculated energy is invariant under rigid deformation 25|, |26] . One of 
the corresponding deformation gradient tensor F z for uniaxial strain in the zigzag direction 
is selected as 

F 

where e z is the stretch ration e in the zigzag direction, e is determined by the Lagrangian 
elastic strain through equation 
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-e 2 + e-r] = 0. 
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The stress-strain relationships of the uniaxial strain in the zigzag direction are 
£1 = C 12 r) z + -C ll2 r,l + -C lll2 r,l + —Cmi2vt, 



£2 = CnVz + 7,CmV z + gCim% + ^C 1U nV z , 
For uniaxial strain in the armchair direction, the strain tensor is, 
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One of the corresponding deformation gradient tensor F a for uniaxial strain in the arm- 
chair direction is 
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where e a is e in the armchair direction. The stress-strain relationships are 



Si = C n r] a + -C 22 2vl + qC 2 222vI + -^C 22222 rf a , 



(10) 



n = C 12Va + i(Cin - C 222 + C ll2 )rf a 

+ t^(Ciiii + 2Cni2 — C 2222 )i] a + —Ci 2222 i] a , 



(11) 



For equibiaxial strain in-plane, r\ a — Vz — 77, the strain tensor is, 
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The corresponding deformation gradient tensor F& for equibiaxial strain in-plane is 
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(13) 



. The stress-strain relationships are 
E* = E 6 
+ 
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All fourteen elastic constants contribute to the expressions for stress-strain response for these 
three deformation states. But the components of SOEC, TOEC, FOEC are over-determined. 
As the second step, we apply the least-square solution to these over- and well- determined 
linear equations. We used least-squares solutions to solve the equations A • C = X by 
computing the elastic constants that minimizes the Euclidean 2-norm ||S — A • C|| 2 . For 
SOEC components, Cn,Ci 2 is obtained by 
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where Ef (Oi) is the coefficient of the first order of strain in E^ (Eqn.9). Similar notation for 
the others. The Young's modulus is E — {C\ x —C\^)jC\\ and Poisson's ration is v — Cu/Cn. 
For TOEC components Cm, Cn 2 cm<iC 222 are obtained by 
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For FOEC components Cim, Cm 2 , Cn 22 and C 2222 are obtained by 
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For FFOEC components Cum, Cum, C n i22, C12222 and C22222 are obtained by 
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III. DENSITY-FUNCTIONAL CALCULATIONS 

The stress-strain relationship of graphene under the desired deformation configurations 
is characterized via ab initio calculations with the density-functional theory (DFT). DFT 
calculations were carried out with the Vienna Ab-initio Simulation Package (VASP) |27h 
30] which is based on the Kohn-Sham Density Functional Theory (KS-DFT) 3l|, |32| with 
the generalized gradient approximations as parameterized by Perdew, Burke and Ernzerhof 
(PBE) for exchange-correlation functions [33]. The electrons explicitly included in the cal- 
culations are the (2s 2 2p 1 ) electrons of boron and (2s 2 2p 3 ) electrons of nitrogen. The core 
electrons (Is 2 ) of boron and nitrogen are replaced by the projector augmented wave (PAW) 



and pseudo-potential approach [34 



35J . A plane- wave cutoff of 520 eV is used in all the 



calculations. The calculations are performed at zero temperature. 

The criterion to stop the relaxation of the electronic degrees of freedom is set by total 
energy change to be smaller than 0.000001 eV. The optimized atomic geometry was achieved 
through minimizing Hellmann-Feynman forces acting on each atom until the maximum forces 
on the ions were smaller than 0.001 eV/A. 

The atomic structures of all the deformed and undeformed configurations are obtained by 
fully relaxing a 24-atom-unit cell where all atoms were placed in one plane. The simulation 
invokes periodic boundary conditions for the two in-plane directions while the displacement 
to out-of-plane direction is forbidden. 

The irreducible Brillouin Zone was sampled with a Gamma-centered 19 x 19 x 1 k- 
mesh. Such large fc-mesh was used to reduce the numerical errors caused by the strain of 
the systems. The initial charge densities were taken as a superposition of atomic charge 
densities. There was a 14 A thick vacuum region to reduce the inter-layer interaction to 
model the single layer system. The results of the calculations are independent of the precise 



8 



value of the out-of-plane thickness, so there is no physical interpretation attached to the 
quantity. 

The VASP simulation calculates the true or Cauchy stress, <r, which for monolayer h-BN 
must be expressed as a 2D force per length with units of N/m by taking the product of the 
Cauchy stress (with units of N/m 2 ) and the super-cell thickness of 14 A. The Cauchy stress 
is related to the second Piola- Kirchhoff(PK2) stress S as 

£ = JF-VtF" 1 ) 11 (19) 

where J is the determinant of the deformation gradient tensor F. 

IV. RESULTS AND ANALYSIS 

We first optimize the equilibrium lattice constant for monolayer h-BN. The total energy 
as a function of lattice spacing is obtained by specifying several lattice constants varying 
around 1.45 Awith full relaxations of all the atoms. A least-squares fit of the energy vs lattice 
constant with a fourth-order polynomial function yields the equilibrium lattice constant, 
aO = 1.4503 , which corre- sponds to the minimum total energy. The result is in good 



agreement with experiments 36[] in h-BN (2.51 A) This most energy favorite structure is set 
as the strain-free structure in this study and the geometry is shown in Fig. [TJ 

When the strains are applied, all the atoms are allowed full freedom of motion within 
plane. A quasi-Newton algorithm is used to relax all atoms into equilibrium positions within 
the deformed unit cell that yields the minimum total energy for the imposed strain state of 
the super cell. 

Both compression and tension are considered here in order to sampling larger elastic 
deformation region. We studied the behavior of the system under the Lagrangian strain 
ranged from -0.1 to 0.3 with a increment of 0.02 in each step for all three cases. There are 
63 ab initio DFT calculations in total. 

The system's energy will increase when strains are applied. Here we define strain energy 
per atom E s = (E tot — E )/n, where E tot is the total energy of the strained system, E is 
he total energy of the strain-free system, n is the number of atoms in the system. Fig. [2] 
shows the E s as a function of strain in uniaxial armchair, uniaxial zigzag and equibiaxial 
deformation. E s responses differently at different strain direction, consistent to the non- 
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FIG. 2. Energy-strain responses for uniaxial strain in armchair and zigzag directions, and equibi- 
axial strains. 
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FIG. 3. Stress-strain responses for uniaxial strain in armchair and zigzag directions, and equibiaxial 
strains. The continuum responses are the least-square fit of the ah initial DFT calculations. 

isotropic structure of the monolayer h-BN. E s are non-symmetrical for compression (77 < 0) 
and tension (77 > 0) for all three cases. This non-symmetry indicates the anhomonicity of 
the monolayer h-BN structures. E s deviated from quadratic relationship with 7/ at strain of 
0.04 in the three tested deformations. 

The stress (second P-K stress) strain (Lagrangian strain) relationship for uniaxial strain 
in armchair and zigzag directions, and equibiaxial strains are shown in Fig. El These stress- 
strain curves reflects the facts of the non-isotropic h-BN structure and anharmonic response 
in compression and tension. For the uniaxial deformation along armchair direction, the 
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TABLE I. Nonzero independent components for the SOEC, TOEC, FOEC and FFOEC tensor 
components, Poisson's ration v and in-plane stiffness Y s of h-BN from DFT calculations. 



SOEC 


TOEC 


FOEC 


FFOEC 


(N/m) 


(N/m) 


(N/m) 


(N/m) 


Cn=293.1 


C m =-2515.8 


diii=18161 


Cnui=-65265 


C 12 =63.76 


Cn2=-428.5 


Clll2^5ooD 


Cmi2=-8454 




C 222 =-2300.6 


Cii22=-2868 


Cni22=-67857 


Y s =279.2 




(^2222 = 12451 


C*12222=-10780 

C 2 2222=-117409 



maximum stress of E" = 5.08, S° = 23.56 (N/m) at r] a = 0.18. For the uniaxial deformation 
along armchair direction, the maximum stress of S^ = 26.26 (N/m) at rj z = 0.26 and 
£2 = 7.82 (N/m) at r] z = 0.30. For the equibiaxial deformation, the maximum stress of 
Sj = T, b 2 = 27.81 (N/m) at 77 = 0.24. 

The elastic constants are the continuum description of the elastic properties. Once we 
know the elastic constants, we can easily apply it for the continuum description. The 
continuum responses are the least-square fit to the stress-stain results from the ab initial 
DFT calculations, as plotted in Fig. El by the equations of ( (ED , (JT]) , (HUD , (IHD , (JHJ) ) . We then 
have 20 values for the fourteen independent elastic constants of h-BN from ab initio DFT 
calculations. The fourteen independent elastic constants of h-BN are finally determined by 
solving equations from Eqn. ( ( 1T5"j) . ( TTB|) . f|T7|) . f[T8~j) ) . The results of these fourteen independent 
elastic constants are grouped in SOEC, TOEC, FOEC and FFOEC and listed in Table 
HI The in-plane Young's modulus Y s = 279.2 (N/m) and Poison's ratio v = 0.2176 were 
obtained from C\\ and Cyi- Our results of C\u C12, Y s and v are comparable with ab initio 
predictioin |37| and tight-binding calculations^ of BN nanotubes. 

The knowledge of these high order elastic constants is very useful to understand the 
anharmonicity. With the high order elastic constants, we can easily study the second-order 
elastic moduli on the pressure p acting in the plane of monolayer of h-BN sheet. Explicitly, 
while the pressure is applied, the second-order elastic moduli are transformed according to 



the relationships 
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C\\ — C\\ — [C\ 
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ClUj-yT-P, 



(20) 
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C-22 — Cu — C222—TP — Pi 
I s 

C\2 = C\2 — Cn2—rp — p. 
is 



(21) 



(22) 



The second-order elastic moduli increase linearly with the applied pressure p within the third- 
order term trucation, as demonstrated in Fig. HJ These equations and plots also indicate 
that the h-BN layer respond to compression or tension along different directions in different 
manners. While pressure presented, the Cu is not symetrical to C22 any more, although 
the difference is relatively small. Olny when p = 0, C\\ = C22 = Cu. This non-isotropy 
behavior could be the outcome of the anharmonicity. 
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FIG. 4. Predicted second-order elastic moduli varies with the pressure p acting in the plane of 
monolayer of h-BN sheet 



V. CONCLUSION 

In summary, we introduced a method to obtain the continuum description of the elastic 
properties of monolayer h-BN through ab initio density functional theory. This thermody- 
namically rigorous continuum description of the elastic response is formulated by expanding 
the elastic strain energy density in a Taylor series in strain truncated after the fifth-order 
term, we obtained a total of fourteen nonzero independent elastic constants for the up to 
tenth-order tensor. We predicted the pressure dependent second-order elastic moduli. This 
continuum formulation is suitable for incorporation into the finite element method. 
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